Periodic magnetorotational dynamo action as a prototype of 
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The nature of dynamo action in shear flows prone to magnetohydrodynamic instabilities is investigated using 
the magnetorotational dynamo in Keplerian shear flow as a prototype problem. Using direct numerical sim- 
ulations and Newton's method, we compute an exact time-periodic magnetorotational dynamo solution to the 
three-dimensional dissipative incompressible magnetohydrodynamic equations with rotation and shear. We dis- 
cuss the physical mechanism behind the cycle and show that it results from a combination of linear and nonlinear 
interactions between a large-scale axisymmetric toroidal magnetic field and non-axisymmetric perturbations am- 
plified by the magnetorotational instability. We demonstrate that this large-scale dynamo mechanism is overall 
intrinsically nonlinear and not reducible to the standard mean-field dynamo formalism. Our results therefore 
provide clear evidence for a generic nonlinear generation mechanism of time-dependent coherent large-scale 
magnetic fields in shear flows and call for new theoretical dynamo models. These findings may off'er important 
clues to understand the transitional and statistical properties of subcritical magnetorotational turbulence. 



PACS numbers: 47.65.Md, 47.20.Ft, 47.27.De, 47.27.Cn 
I. INTRODUCTION 

The generation of coherent, system-scale magnetic 
fields in flows of electrically conducting fluids is a long- 
standing magnetohydrodynamics (MHD) problem which 
has so far mostly been analyzed in terms of linear, kine- 
matic mean-field dynamo action |1|. Mean-field eflTects 
have in particular long been invoked to explain the origin 
of magnetic cycles in MHD rotating shear flows 12) [31 . 
However, there is currently no mathematical theory and 
only little physical understanding of mean-field dynamos 
in parameter regimes (kinetic and magnetic Reynolds num- 
bers) typical of laboratory dynamo experiments and natural 
MHD shear flows. Overall, the physical nature of dynamo 
action in such flows remains a mostly open question, with 
critical implications for geophysics and astrophysics. 

Studying the dynamics of shear flows prone to the de- 
velopment of various local three-dimensional (3D) MHD 
instabilities is one of the most promising avenues of re- 
search on this problem. Numerical simulations of this 
class of flows have explicitly demonstrated their potential 
for coherent dynamo action ll4UT3]l, and it has occasion- 
ally been pointed out that their dynamics diflfers markedly 
from that of kinematic mean-field dynamos, as magnetic 
field amplification does not proceed exponentially in time 
and requires finite- amplitude magnetic perturbations cou- 
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pling dynamically to fluid motions. Identifying the physi- 
cal mechanisms underlying this behaviour is of prime im- 
portance to improve our understanding of dynamo action. 

It has recently been realized that these instability-driven 
Q [H or subcritical dynamos in shear flows (TOl [141 [El 
bear many similarities to the hydrodynamic transition to 
turbulence of shear flows [16-20], a fundamentally non- 
linear process whose dynamics involves a variety of hy- 
drodynamic coherent structures such as equilibria, travel- 
ling waves 1 21 - 24] or limit cycles | 25 -28|. This naturally 
raises the question of the existence and dynamical rele- 
vance of coherent structures in MHD shear flows. Most 
importantly for dynamo theory, could time-dependent 
large-scale dynamo action in such flows be explained in 
terms of 3D nonlinear dynamo cycles ? Besides, if nonlin- 
ear dynamo cycles exist for this kind of systems, what can 
be learned from them regarding how turbulence originates 
and develops in MHD shear flows ? 

The problem of magnetorotational (MRI) dynamo ac- 
tion in Keplerian flow (often referred to as "zero net flux 
MRI"), encountered in the context of astrophysical accre- 
tion disks 1 29 1, is particularly interesting to address these 
questions, as direct numerical simulations display pseudo- 
cyclic magnetic dynamics |4l[rT][T2]| and indicate that tran- 
sition to sustained MHD turbulence is intrinsically 3D and 
nonlinear |5 1. Besides, numerical MRI dynamo turbulence 
seems to have finite lifetime and to be structured around 
a chaotic saddle |[T5ll . much like hydrodynamic turbulence 
in pipe flow | TS] [191. The search for coherent MRI dy- 
namo structures such as nonlinear cycles is still in its in- 
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fancy, though. Progress on these matters is not only de- 
sirable from a general dynamo perspective, it would also 
shed light on the transitional properties of MRI turbulence 
|[30l and on its dynamical properties (saturation, transport) 
in fully developed regimes. 

The only exact MRI dynamo structure known as yet is 
a nonlinear equilibrium in Keplerian MHD plane Couette 
flow with walls ifTOl . and it only exists for a limited range 
of parameters. In this work, we report the first accurate 
calculation of a 3D nonlinear MRI dynamo limit cycle 
at moderate (transitional) kinetic and magnetic Reynolds 
numbers using numerical techniques similar to those used 
for computing nonlinear hydrodynamic cycles in plane 
Couette flow |27|. The study of this nonlinear MRI dy- 
namo cycle enables us to investigate in detail the essential 
physical mechanisms underlying sustained time-dependent 
MRI dynamo action. In particular, we show that the cycle 
dynamics is not amenable to a standard mean-field dynamo 
description, thereby providing clear and detailed evidence 
for a fully nonlinear mechanism of coherent magnetic field 
generation in shear flows prone to MHD instabilities. 

The paper is organized as follows. In Sect. |llj we in- 
troduce the theoretical framework and numerical methods 
used in the study. In Sect. [IIl| we discuss the physical 



principles of the MRI dynamo and describe our strategy 
to excite large-scale, coherent recurrent dynamics in direct 



numerical simulations. Section IV is devoted to the pre- 
sentation and detailed analysis of the main result of the pa- 
per, namely the discovery of a nonlinear MRI dynamo cy- 
cle computed thanks to a Newton algorithm. A discussion 
(Sect. |V]) of the results and of their possible implications 
for dynamo theory, astrophysics and research on hydrody- 
namic shear flows concludes the paper. 




X (radial) 

FIG. 1 : Geometry of the shearing sheet. For a Keplerian flow, the 
vorticity of Us is anti-aligned with the rotation vector. 



sipative MHD equations in an unstratified shearing sheet: 
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V-M = 0, V-B = 0. 



(3) 



The kinetic and magnetic Reynolds numbers are defined 
according to Re = SL^/v and Rm = SL^/rj, where v and 
T] are the constant kinematic viscosity and magnetic diffu- 
sivity and L is a typical scale of the spatial domain consid- 
ered, n is the total pressure (including magnetic pressure) 
and B is expressed as an equivalent Alfven velocity. In 
the following, velocity and magnetic field amplitudes are 
measured with respect to S L, while times are measured 
with respect to the inverse of the shearing rate S~^. 



II. THEORETICAL AND NUMERICAL FRAMEWORK 



B. Strategy for capturing cycles 



A. The shearing sheet 

We use the shearing sheet approach to diff'erentially ro- 
tating flows 1 31 1, whereby a cylindrically symmetric dif- 
ferential rotation profile is approximated locally by a linear 
shear flow Us = -Sx ey and a uniform rotation rate O = 
^e^. For a Keplerian flow, Q = (2/3)5'. Here, (x,y,z) are 
respectively the shearwise, streamwise and spanwise di- 
rections (radial, azimuthal and vertical in accretion disks). 
The geometry of the shearing sheet is represented in Fig.[T] 
To comply with dynamo terminology, we refer to the (x, z) 
projection of vector fields as their poloidal component and 
to their y projection as their toroidal component. "Axisym- 
metric" fields have no y dependence. For readers familiar 
with hydrodynamic plane Couette flow, non-axisymmetric 
perturbations in our problem correspond to "streamwise- 
dependent" perturbations in plane Couette flow. 

We consider incompressible velocity perturbations u 
and magnetic field B whose evolutions obey the 3D dis- 



Eqs. ([TJ-|3 often implemented numerically using 
the cartesian shearing box model described in Sect. |IIC 
below. In the regime of Keplerian diff'erential rotation, 
MRI dynamo action has been found in a large number 
of independent direct numerical simulations of this kind 
El El [TT] [53^1, sometimes in the form of pseudo-cyclic 
dynamics. Our strategy to capture a nonlinear MRI dy- 
namo cycle in this regime follows that used in Ref. 1271 
to compute nonlinear cycles in hydrodynamic plane Cou- 
ette flow. First, we tried to excite pseudo-cyclic dynam- 
ics in direct shearing box numerical simulations (DNS) by 
devising initial conditions inspired by our partial knowl- 
edge of the underlying physics. Once this was achieved, 
we attempted to converge to a cycle using a Newton solver 
seeded with a well-chosen DNS snapshot and an estimate 
for the cycle period. The first part of this program is de- 
scribed in Sect.lnil Convergence to a cycle using Newton's 



method is presented in Sect. IV The numerical methods 
for achieving this result are explained below. 
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C. Numerical method for direct numerical simulations 

The numerics presented in this paper are done in the in- 
compressible cartesian shearing box framework, which as- 
sumes simple spatial periodicity in the y and z directions 
and shear-periodicity in the x direction. The latter amounts 
to assuming that a linear shear flow is constantly imposed 
and that all quantities are periodic in a sheared Lagragian 
frame. This approximation is justified in the context of 
centrifugally supported diff'erentially rotating flows, such 
as as trophy sical accretion flows. Simulations of this kind 
are also sometimes referred to as homogeneous shear flow 
turbulence simulations (see e.g. Ref. |l33l,^34l). 

Time integrations (direct numerical simulations) are car- 
ried out in a shearing box of size (Z^x? ^y-> ^z) using the 
SNOOPY code |35|. The numerical time integration 
scheme used is a standard explicit third-order Runge-Kutta 
algorithm. The code relies on a spectral implementation 
of the shearing box model similar to that described in 
Ref. 1361 . A discrete spectral basis of "shearing waves" 
(or Orr-Kelvin waves, see Refs. (371 [38l) with constant ky 
and kz wavenumbers and constant shearwise Lagrangian 
wavenumber kf^ is used to represent the various fields 
in the sheared Lagragian frame. The shearing of non- 
axisymmetric perturbations in this model is described us- 
ing time-dependent Eulerian shearwise wavenumbers, 

k^^kf +kySt . (4) 

This equation for the radial wavenumber provides an ex- 
act description of the evolution of non-axisymmetric waves 
with initially leading polarization {kf"^ ky < 0) into trailing 
waves (kxif) ky > 0, corresponding to a trailing spiral in 
cylindrical geometry) under the action of shear. 

An important comment is in order at this stage. If we 
were to simply consider the evolution of a given initial 
set of such waves, we would only be able to observe dy- 
namics that decay on long times. Indeed, the fate of all 
shearing waves is to evolve into strongly trailing structures 
with ever smaller scale in x, and such structures are ex- 
tremely efficiently dissipated by viscous and resistive pro- 
cesses (see e.g. Ref. ||39l|40l). In order to accomodate for 
possible physical interactions leading to long-lived nonlin- 
ear dynamics in numerical shearing box simulations, a pro- 
cedure must therefore be used that leaves open the possi- 
bility of physical generation and dynamical evolution of 
new leading non-axisymmetric structures in the course of 
the simulations. The solution to this problem is to regu- 
larly "remap" the basis of shearing waves used to describe 
the various fields. At regular time intervals during the 
simulations, the energy content of strongly trailing shear- 
ing waves is set to zero, the corresponding basis vector is 
pruned and replaced by a new shearing wave basis vector 
with strongly leading wavenumber (see Ref. |41 1, Chap. 5, 
Sect. 4). If the simulation is well resolved spatially (as is 
the case for the results presented in this manuscript), the 
energy contained into strongly trailing waves when they 
are pruned should be negligible (as a result of their en- 



hanced dissipation), and the remap procedure should not 
therefore artificially aff'ect the dynamical evolution of the 
system in any significant way. A way to check this is to 
compare the energy lost by this procedure with the en- 
ergy dissipated by viscous/resistive diff'usion. We always 
find that artificial energy losses are negligible in SNOOPY 
for spatially well-resolved simulations |42|. We also point 
out that the remap procedure does not by itself inject en- 
ergy into new leading waves but merely provides room for 
them in the wavenumber grid. Finally, we emphasize that 
all nonlinearities of the shearing sheet MHD Eqs. ([T])-([3]), 
including nonlinear interactions between all the shearing 
waves present at a given spatial resolution, are retained in 
our numerical model. A standard pseudo- spectral method 
with dealiasing is used to compute all nonlinear terms at 
each time step. 

D. Newton's method for computing nonlinear cycles 

Newton's method is a standard tool for computing non- 
linear coherent structures such as saddle points, travelling 
waves, or nonlinear cycles in high-dimensional dynami- 
cal systems. In recent years, the method has been applied 
successfully to the three-dimensional Navier-Stokes equa- 
tions for various wall-bounded shear flows ll2TH24l[27ll28l 
and to the MHD equations in Keplerian plane Couette flow 
ifTOl . For the purpose of this study, we developed a new 
Newton solver called PEANUTS. The solver makes use of 
the PETSc toolkit |43 1 and is based on an efficient matrix- 
free Newton-Krylov algorithm particularly well adapted 
to calculations for high-dimensional dynamical systems 
such as those resulting from the discretization of the three- 
dimensional partial diff'erential equations of fluid dynam- 
ics. It can be used to compute nonlinear equilibria, travel- 
ling waves and limit cycles for a variety of partial difl'eren- 
tial equations. For a nonlinear cycle search, the code mini- 
mizes ||Z(r)-Z(0)||2/||Z(0)||2, where X{t) is a state vector 
containing all independent field components at time t, and 
r is a guess for the period 127 1 . An eigenvalue solver based 
on the SLEPc toolkit 1441 was implemented to compute the 
stability of nonlinear states. The code was tested against 
solutions to the Kuramoto-Sivashinsky equation |45 | be- 
fore being implemented for the 3D MHD equations in the 
shearing box, using SNOOPY as time integrator. 

III. EXCITATION OF RECURRENT DYNAMICS 

In this Section, we describe our strategy to approach a 
nonlinear MRI dynamo cycle using DNS of the 3D MHD 
equations in the shearing box. We first discuss in detail 
what is the "minimal" set of initial conditions required to 
excite a long-lived MRI dynamo in direct numerical sim- 
ulations. We then explain how smooth-enough pseudo- 
cyclic dynamics can be excited at moderate Re and Rm 
with this kind of initial conditions by varying the aspect 
ratio of the simulations and by restricting the dynamics to 
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an invariant subspace associated with a natural symmetry 
of the original equations. 



A. Devising a good initial guess for a Newton search 

Previous work has demonstrated that instability-driven 
dynamo action requires a dynamical interplay between 
a "large-scale" axisymmetric, instability-supporting mag- 
netic field and perturbations unstable to non-axisymmetric 
MHD instabilities, whose amplification to nonlinear levels 
may generate an electromotive force (EMF) with the abil- 
ity to sustain the large-scale field l7 tiTTl[nifT4ll46ll . The 
basic processes thought to be responsible for MRI dynamo 
action are described below and in Fig. |2] Let us focus 
on the time-evolution of the axisymmetric magnetic field 
B(x, z, t), where the overbar denotes an average over y, and 
JqJq'' B(x, z, t) dxdz = for the MRI dynamo ("zero net- 
flux") problem. The induction equation for B reads 



dB 

~dt 



-S B^ey-\-V xE -\-T]AB . 



(5) 



The first term on the r.h.s. describes the stretching of the 
axisymmetric poloidal field into an axisymmetric toroidal 
field (the so-called D eff'ect in dynamo theory), the second 
term is a nonlinear induction term involving the axisym- 
metric projection E = ux B of the electromotive force 
resulting from the nonlinear coupling of velocity and mag- 
netic perturbations, and the third term is the magnetic dif- 
fusion term. In the following, we will be particularly in- 
terested in the time evolution of the fundamental Fourier 
mode inzof B, defined as 



Bo(z, t) = Bo(t)cos (k,oz) 



(6) 



with = 2n/Lz, as this mode is always and by a large 
amount the dominant contribution to the total axisymmet- 
ric field for the type of dynamics excited by the class of 
symmetric initial conditions described below (a large-scale 
field with an arbitrary phase in z can of course be excited 
if the simulation is initialized with non- symmetric pertur- 
bations or if symmetry-breaking instabilities are allowed 
to develop during the simulations). Following Eq. ([5]), we 
also introduce the nonlinear EMF acting on 



Eo(z, t) = Eo(t)sm(k,oz) . 
The time-evolution of Bo(t) is given by 



(7) 



dBo 

dt 



-S BoM ey + ^,0 e, x Eo(t) - r]k%Bo(t) . (8) 



The physical interpretation of the various terms on the r.h.s. 
of this equation is the same as for Eq. ([5]). Our goal to ob- 
tain sustained time-dependent MRI dynamo action in di- 
rect numerical simulations was to excite a dynamo loop 
involving the first two terms on the r.h.s. of Eq. ([8]), as de- 
picted in Fig. [2] To do so, we first attempted to start from 
a very simple class of initial conditions combining 
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FIG. 2: (Color online) Suggested physical mechanism of the 
MRI dynamo. Full arrows: main dynamo loop. Dashed ar- 
rows: nonlinear regeneration of MRI-unstable fluctuations. The 
various colors are used to identify the active modes taking part 
in the cyclic MRI dynamo described in Fig. [5] red (top left 
box) and blue (top right box) denote axisymmetric field com- 
ponents, diflferent colors in the bottom boxes denote successive 
non-axisymmetric MRI waves. 



1. a "large-scale", axisymmetric poloidal magnetic 
field ^ojc (blue box in Fig. [2]) whose stretching by the 
shear ("Q eff'ect", full line blue arrow in Fig.|2]) pro- 
duces an MRI-unstable axisymmetric toroidal field 
Boy (red box in Fig. [2]), 

2. non-axisymmetric perturbations subject to joint am- 
plification (full line red arrow in Fig. [2]) by a tran- 
sient toroidal MRI |46 -49 | of Boy and by a kine- 
matic Orr mechanism (swing amplification of non- 
axisymmetric waves, see Refs. |[381|Ml|2])- These 
perturbations were chosen in the form of real ran- 
dom shearing wave packets with a full spectrum of 
kz and a single "horizontal" wavenumber pair (kf\ 
ky) = ±(-2n/L^,2n/Ly) with leading (kfhy < 0) 
polarization (diff'erent colors in the bottom left box 
of Fig. [2] represent successive individual leading 
waves). In the course of their evolution, such per- 
turbations may generate a nonlinear electromotive 
feedback ^0 with the ability to sustain (full line 
green arrows in Fig. [2]), thereby closing the main dy- 
namo loop. 

It turns out that, independently of the initial amplitudes 
of each of these perturbations, this restricted class of initial 
conditions can only trigger transient, short-lived dynamics. 



As explained in Sect. II C non-axisymmetric perturbations 
with a single initial horizontal wavenumber (kf\ ky) can 
only be amplified for a few shearing times before they get 
sheared into a strongly trailing (^^(0 ky » 0, bottom right 
box with multiple colors in Fig. [2]), rapidly decaying struc- 
ture 1461 . Besides, their nonlinear self-interaction cannot 
give rise to new non-axisymmetric leading waves, making 
it impossible to sustain against ohmic diff'usion on long 
times. Hence, long-lived dynamics can only be excited if 
a distinct physical mechanism operates that generates new 
leading, transiently MRI-unstable perturbations. Explor- 
ing this issue, we then found that much longer-lived dy- 
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namics is obtained as soon as 

3. an x-dependent axisymmetric modulation 

^mod(-^, t) = 5mod(0 COS {271x1 L^) COS (^^o^) (9) 

of B is initially added on top of i-c 
B{X, Z, ^ = 0) = Bq{z, ^ = 0) + ^mod(^, z, ^ = 0) . (10) 

This simple numerical observation indicated that the phys- 
ical mechanism by which new leading shearing waves are 
generated requires that B be modulated along the x direc- 
tion. The physical explanation for this behaviour is rather 
subtle, though : such a modulation of B confines MRI- 
unstable perturbations in x and allows for reflections of 
non-axisymmetric waves, somehow taking on the role of 
walls in a wall-bounded shear flow (non-axisymmetric in- 
stabilities in wall-bounded shear flows take on the form 
of global standing modes, see e.g. Ref. |17|). In more 
mathematical terms, the nonlinear triad interaction (shown 
by dashed arrows in Fig. [2]) of the "confining" (say with 
kx = -27i/Lx,ky = 0) axisymmetric mode ^mod (red 
dashed arrow) with a trailing shearing wave (say kx(t) = 
71 /Lx, ky = 271 1 Ly, green dashed arrow) can seed a new 
leading (^^(0 = -^/^x, ky = 2nlLy) wave (light blue 
dashed arrow). This type of mechanism, which has long 
been suspected to be at work in non-rotating hydrodynamic 
shear flow turbulence (see for instance the discussion in 
Ref. 1521 ) and has more recently been invoked in the con- 
text of hydrodynamic stability of accretion disks | 53 1, is in 
fact essential to any sustained non-axisymmetric dynamics 
in the shearing box. We did check that the seeding of new 
leading waves in the simulations was not related to our im- 
plementation of the remap procedure (see Sect. |IIC| and 
Ref. (36J), but had a genuine physical origin. In particular, 
the fact that long-lived dynamics and new leading waves 
can only be excited in the simulations if a ^mod compo- 
nent is added to the initial condition demonstrates that our 
numerical method does not artificially inject energy into 
leading waves. 

To summarize this paragraph, using various relative 
combinations of axisymmetric and non-axisymmetric ini- 
tial conditions composed of ^o, ^mod and a random lead- 
ing shearing wave packet, we found it possible to obtain 
long-lived MRI dynamo action in direct shearing box sim- 
ulations for diflferent regimes. Based on these experiments, 
we claim that the essence of the driving mechanism of the 
MRI dynamo can be fully explained in terms of the few 
generic physical mechanisms described above and in Fig- 
ure [2] As will be shown below, this claim is well supported 
by the targeted numerical experiment presented in Sect.|lV] 

B. Simplifying the dynamics 

Approaching nonlinear cycles by DNS finally required 
us to find regions of parameter space in which only a small 
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FIG. 3: (Color online) Projection in the Bq^ - Bqj plane of a 
DNS trajectory approaching a nonlinear cycle for {Ly^,Ly,L^) - 
(0.1,20,2) L and Re = 70, Rm = 360. The full blue circle 
whose coordinates are (0.046, 0) marks the position of the sys- 
tem at ^ = and the red one located at (-0.046, 0.997) marks 
the position at ^ = 457.1 S~^. Typical amplitudes of the various 
components of the initial conditions used to obtain this kind of 
trajectory are given in Sect. |III C] 

number of these structures are present. The dynamics in 
regimes (Re, Rm of a few thousands) typical of simula- 
tions displaying pseudo-cyclic dynamics (H [TT| [121 be- 
ing complex and probably involving a lot of diff'erent co- 
herent structures, we restricted our investigations to Re, 
Rm of a few hundreds. In such regimes, however, shear- 
ing waves are quickly damped after they turn trailing, un- 
less kx(t) changes on a timescale much longer than S~^. 
Exciting long-lived dynamics in this Re and Rm regime 
therefore further required setting Ly » Lx (\ky\ «c |^^^^| 
in Eq. (|4])). Starting from various non- symmetric initial 
conditions, we then spotted that nonlinear states approach- 
ing a symmetry | 54 | of the shearing box MHD equa- 
tions were regularly excited in DNS. This symmetry, de- 
scribed in the Appendix, allows for a large-scale axisym- 
metric magnetic field with the symmetry of ^o- In order to 
isolate these structures more easily, we then enforced nu- 
merically that the dynamics take place in the correspond- 
ing invariant subspace. This strategy was sufficient to ex- 
cite recurrent dynamics in a large aspect ratio shearing box 
with (Lx,Ly,L,) = (0.7, 20, 2) L and Re = 70, Rm = 360. 
A projection in the Bqx - B^y plane of a DNS trajectory 
approaching a nonlinear cycle in this regime is depicted in 
Fig. |3] The trajectory is seen to approach a periodic orbit 
after a few tens of shearing times and then to stay close to 
it for several hundred shearing times. 

C. Sensitive dependence on initial conditions 

Before we close this Section, we find it necessary to em- 
phasize that the dynamical system at hand has a very high 
sensitivity on initial conditions. This feature seems to be 
a very common property of shear flow turbulence (see for 
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instance Ref . fTS] for the case of hydrodynamic turbulence 
in pipe flow) and was explicitly demonstrated for the MRI 
dynamo problem in Ref. ITSll . In the course of this work, 
we clearly observed that the set of initial conditions lead- 
ing to long-lived dynamics for the system at hand is fractal. 
For this reason, any statement of a set of precise numerical 
values for the amplitudes of the various types of perturba- 
tions involved in the design of the initial conditions leading 
to pseudo-cyclic dynamics is rather pointless. The prob- 
lem is that difl'erent numerical codes, initialized exactly 
in the same way, will almost certainly diverge after a few 
shearing times because their algorithms will generate dif- 
ferent numerical "noise" (low amplitude numerical errors) 
at each time step, subsequently leading to a rapid diver- 
gence (between difl'erent codes and also possibly difl'erent 
computer architectures) of phase- space trajectories such as 
that shown in Fig. |3] The initial amplitudes and relative 
mixtures of modes required to approach nonlinear cycles 
are therefore in the end specific to the numerical meth- 
ods used in the code (this does not imply that the cyclic 
nonlinear solutions are not themselves robust, as discussed 
in the next Section). The only relevant helpful informa- 
tion for a reader eager to reproduce the results presented in 
Sect. [IV| and trajectories similar to that shown in Fig. [3] is 
a set of approximate values for the amplitudes of the per- 
turbations entering our class of initial conditions, around 
which we found it possible to excite recurrent dynamics 
for (L^, Ly, L,) = (0.7, 20, 2) L and Re = 70 and Rm = 360: 
Bqx - 0.0465'L, ^mody - 0.115'L, and non-axisymmetric 
shearing wave packets (with random relative amplitudes 
for diff'erent k^, as described earlier) with comparable ve- 
locity and magnetic amplitudes of the order 03SL. 



IV. A NONLINEAR MRI DYNAMO CYCLE 



Following the strategy defined in Sect. |IIB[ we then at- 
tempted to capture precisely the nonlinear cycle underlying 
the recurrent dynamics spotted in Fig.[3]using the Newton- 
Krylov solver PEANUTS. 



A. Convergence with Newton's method 

Because of shearing-periodic boundary conditions, only 
cycles whose period is a multiple of Ly/(SLx) are al- 
lowed in shearing boxes. The periodicity of the recur- 
rent dynamics spotted in the DNS being very close to 
T = 2Ly/(SLx) ^ 57.143 this period was imposed on 
the Newton solver. Initializing the solver with a snapshot 
of a DNS displaying dynamical recurrences (the red full 
circle point of Fig. [3]), we obtained convergence to a cycle 
after 0(10) iterations with a relative error of 10"^. Each 
Newton step requires 0(10) Krylov iterations to solve the 
Jacobian system with a relative error smaller than 10"^. 
As each Krylov iteration requires running a DNS over a 
cycle period, (9(100) numerical integrations are needed to 



"capture" this cycle with this iterative method. The dy- 
namics of this nonlinear MRI dynamo cycle is illustrated 
in Fig. [4] The results presented here are for a resolution 
(Nx,Ny,Nz) = (48,24,72) with 2/3 dealiasing (the same 
resolution used to approach the cycle by DNS). Half that 
resolution already ensures convergence to very good ac- 
curacy. Analyzing the energy spectrum of the cycle, we 
only found diff'erences of a few percent between the full- 
resolution and half-resolution runs. 



In respect of the observations made in Sect. IIIC 



we 

emphasize that the cycle calculated by the means of New- 
ton's method is not a spurious numerical feature. The 
most important requirement for convergence, as is stan- 
dard with Newton's algorithm, is that the DNS snapshot 
used as a starting guess for the algorithm be in close- 
enough vicinity of the cycle. Provided that this require- 
ment is satisfied, convergence to the cycle is robust and 
can be obtained by using (as a starting guess for the New- 
ton solver) various adequate DNS snapshots resulting from 
integrations of various sets of initial conditions, for var- 
ious Courant-Friedrichs-Lewy conditions for the numeri- 
cal time-integrator, using either fixed or adaptative time- 
stepping. 

Finally, it is important to stress that the cycle at hand 
is a genuine solution to the full discretized nonlinear 3D 
MHD equations in the shearing box, not just a solution to 
a low-dimensional, reduced dynamical model of the prob- 
lem. The imposition of symmetries in the numerical reso- 
lution only helps to target nonlinear cycles with a given 
natural symmetry (i.e. a symmetry allowed by the full 
MHD equations in the geometry studied) more easily. In 
fact, once convergence is achieved with the Newton solver, 
the cycle can be integrated for several periods in a standard 
DNS without enforcing any symmetry. 



B. Description of the cycle 

The physics of the cycle is best understood by looking 
at Fig.kl A large-scale axisymmetric toroidal field By with 
z-dependent polarity dominates at ^ = 0. The MRI me- 
diated by that field progressively amplifies weak, leading 
non-axisymmetric perturbations. The eff'ect of the instabil- 
ity is to separate "radially" (in x) fluid particles initially at- 
tached to individual, almost frozen-in field lines |[29ll , lead- 
ing to their non-axisymmetric distortion (r/8 to T/4). In 
addition, for the cycle at hand, linear MRI velocity per- 
turbations drag fields lines with opposite polarities in op- 
posite directions. As they get sheared, MRI-amplified ve- 
locity perturbations take the form of a non-axisymmetric 
pattern of counter-rotating poloidal flow cells. The z com- 
ponent of the flow advects distorted magnetic field lines 
with opposite polarities in opposite z directions, eff'ec- 
tively reversing regions of positive and negative polari- 
ties (T/4). The X component of the cellular flow even- 
tually advects field lines back to their original x location 
(3r/8), producing at T/2 an axisymmetric magnetic field 
opposite to the original one. Hence, the reversal is the 
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FIG. 4: (Color online) Volume renderings of By and isosurfaces of B - 0.9 colored by By every r/8 (positive By in red/light gray, 
negative By in blue- violet/dark gray). The coordinates of the bottom left comer of each box are x = 0, = 0, z = 0, where the y' 
coordinate is defined as - y - (Lx/2)St (see Appendix). The arrows field in the y' = plane traces non-axisymmetric MRI velocity 
perturbations (velocity perturbations in the y - Ly/2 plane have opposite sign), whose eff'ect is to distort and advect magnetic field 
lines of opposite polarities in opposite directions. 



pure outcome of the nonlinear evolution (self-interactions) 
of non-axisymmetric MRI perturbations. As a result of 
the confinement of MRI perturbations by the shearwise- 
modulated MRI- supporting toroidal field, new leading per- 
turbations seeds are generated during the first half of the 
cycle. Their subsequent amplification and nonlinear evo- 
lution during the second half of the cycle result in a new 
field reversal after another T/2, back to the initial state. 

This qualitative physical scenario is fully supported by 
a quantitative analysis of the time-evolution (Fig. [5]). The 
linear MRI of the z-dependent axisymmetric toroidal field 
Boy (Fig. [5^ in red) amplifies successive non-axisymmetric 
shearing MRI waves (Fig. |5]l, rainbow colors) with ky = 
±271 /Ly and a kz spectrum. The nonlinear self-interaction 
of a single wave packet translates into an axisymmetric 
EMF ^0 whose x and y components are clearly responsible 
for the reversals of (Fig- 1^})). Interestingly, a calcula- 
tion of the various terms on the r.h.s. of the projection of 
Eq. ([8]) on Cy throughout the cycle shows that the Q eflTect 
is not the dominant term in this equation for this particu- 
lar cycle and geometry (Fig. pp). Hence, Bq^ and B^y are 
almost in antiphase, whereas they are almost in quadrature 
for the pseudo-cycle described in Ref. 1 1 1 1. Fig.[5]i clearly 
demonstrates the periodic seeding of new leading pertur- 
bations (yellow, green, light blue etc. represent successive 
shearing waves, as in Fig.|2]). The evolution of the ampli- 
tude of ^modj. defined in Eq. ([9]), is depicted in Fig.pt. 



Bmody is vcry small compared to B^y, which demonstrates 
that a very weak confinement of the MRI is actually sufift- 
cient to generate new leading waves. 

As seen in Fig. |4j the reason why such a cycle can be 
computed at fairly low resolutions (in particular in the y 
direction) is that it is a large-scale, time-dependent coher- 
ent structure whose non-axisymmetric dynamics (and the 
corresponding axisymmetric dynamo feedback ^0 that it 
induces) is largely dominated by MRI shearing waves sup- 
ported by the fundamental Fourier mode ky = 27i/Ly of the 
shearing box in the y direction. We also note that all the 
physical processes involved are not actually specific to the 
shearing box, and may therefore also be present in cylindri- 
cal geometry, with non-axisymmetric perturbations taking 
on the form of sheared spiral waves. 

Finally, remark that the energy required to sustain the 
three-dimensional cyclic dynamics against resistive and 
viscous dissipation is extracted from the shear (which is 
the only available energy source of the system) thanks to 
the toroidal MRI of non-axisymmetric shearing waves. 



C. Stability 

Using the Floquet eigenvalue solver, we found that the 
cycle has a single unstable eigenmode, with a Floquet 
eigenvalue A ^ 32.33 corresponding to a positive growth 
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FIG. 5: (Color online) From top to bottom: evolution over two 
cycles of a) Bqx (blue/dark gray) and ^03; (red/light gray) ; b) -£^03; 
(blue/dark gray, source term for Bq^) and Eq^ (red / light gray, 
source term for ^03;) ; c) amplitudes of the three terms on the r.h.s. 
of the projection of the axisymmetric induction Eq. ^ on ey (dot- 
ted blue: Q effect, full line red: nonlinear EMF, dashed magenta: 
magnetic diffusion) ; d) total energy of successive MRI shear- 
ing waves with ky = ±2nlLy (rainbow colors, dashed/full line: 
leading/trailing phase. Each color represents a different shearing 
wave) ; e), amplitude of the axisymmetric magnetic field modu- 
lation Bmodj- 



rate A = \nA/T ^ 0.061 S (at half resolution, A = 22.63, 
A ^ 0.055 5), and the same ^1 symmetry. Hence, any 
small perturbation to the cycle, independently of its initial 
amplitude, ultimately kicks the system out of its initially 
periodic trajectory in phase- space. This exponential in- 
stability of the cycle was observed in the simulations and 
leads to complete escape from the neighbourhood of the 
cyclic solution after a few hundred shearing times. Unsta- 
ble periodic orbits such as this one are a typical feature of 
chaotic dynamical systems and play an important role in 
their dynamics (see discussion in Sect.jV]). 




-1.0 -0.5 0.0 0.5 1.0 




-0.06 -0.04 -0.02 0.00 0.02 0.04 0.06 



FIG. 6: (Color online) Projection of the periodic MRI dynamo 
orbit in the (blue/dark gray) | Ey (red/light gray) vs. \ By 
planes. 



D. Investigating the nature of dynamo action 

One may finally wonder if the nonlinear couplings lead- 
ing to this dynamo can be described by standard mean- 
field theory, or dynamos are ruled out: both net 
kinetic and current helicities are negligible in the DNS. 
For our cycle, the axisymmetric field and EMF are largely 
dominated by their projection and ^0 on cos(27rz/Lz) 
and sin(2nz/Lz) planforms (see Eqs. ([6])-(|7])), so one may 
be tempted to interpret the dynamo feedback in terms of 
a non-diagonal "turbulent resistivity" fj tensor 1551 that 
would couple the amplitudes of the components of and 
^0 finearly. However, not only is the instantaneous rela- 
tionship between these two quantities unambiguously non- 
linear (Fig. |6]), it also looks rather difficult to model ana- 
lytically based on simple physical arguments. Hence, this 
periodic MRI dynamo does not reduce to a simple standard 
mean-field dynamo, a conclusion that also seems to apply 
to the magnetic buoyancy-driven dynamo 1.8. .56.1. 

V. DISCUSSION AND CONCLUSIONS 

We have presented the first accurate numerical determi- 
nation of a cyclic nonlinear MRI dynamo solution to the 
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original 3D dissipative incompressible MHD equations in 
the regime of Keplerian differential rotation. Preliminary 
investigations indicate that its existence is not limited to a 
narrow range of Re and Rm, unlike the plane Couette flow 
stationary solution reported in Ref. 1 10 |, and may notably 
extend down to low magnetic Prandtl number regimes in 
which recent studies have found it difficult to obtain sus- 
tained MRI dynamo turbulence 1301 . A detailed parametric 
study of the cycle and an assessment of its importance for 
the MRI dynamo transition problem is currently underway 
and will be presented in a separate paper. 

Overall, the discovery of this MRI dynamo limit cycle 
significantly extends and consolidates earlier findings 1 10] 
and claims 1141 03 that dynamo action and MHD turbu- 
lence in shear flows prone to local MHD instabilities has a 
similar nature to hydrodynamic turbulence in shear flows. 
The results notably provide clear evidence that magnetic 
field generation and sustenance at moderate Re, Rm in 
a numerical set-up comparable to that of most MRI dy- 
namo simulations so far results from a nonlinear mecha- 
nism of interaction between axisymmetric fields and non- 
axisymmetric instability modes. The detailed description 
of the process may guide future studies of dynamos medi- 
ated by MHD instabilities other than the MRI |7 - 9 , 1 3 , 53 . 

Our findings also demonstrate that new theoretical mod- 
els of large-scale dynamo action are required. Progress 
on this problem may be possible thanks to periodic or- 
bit theory |58|, which offers formal mathematical connec- 
tions between the statistical properties of dynamical sys- 
tems with a small number of active degrees of freedom 
and their unstable cycles. One should therefore attempt 
to identify new MRI dynamo cycles, analyze their depen- 
dence on various parameters and assess their relative dy- 
namical importance using stability analysis. It is almost 



certain that, similarly to hydrodynamic shear ffows f27l, 
many such structures exist in the phase space of MHD- 
unstable shear ffows at moderate and large Re and Rm. 

This approach may also prove useful for the astrophysics 
problem of turbulent transport of angular momentum in ac- 
cretion disks 1291 . The fact that the MRI dynamo cycle de- 
scribed above generates an average turbulent MHD stress 
a ~ 0.025 directly comparable to that obtained in 

direct numerical simulations of similar conffgurations 1301 
constitutes a notable preliminary ffnding in this respect. 

We ffnally point out that the MHD ffndings reported in 
this study, which have partly been inspired by recent ad- 
vances in the ffeld of shear ffow hydrodynamics, may in re- 
turn be helpful to make progress on understanding hydro- 
dynamic shear ffow turbulence. The leading wave regen- 
eration mechanism that was clearly identiffed in the shear- 
ing box in this work offers an interesting connexion be- 
tween three-dimensional turbulence in wall-bounded and 
homogeneous shear ffows, which may notably help to un- 
derstand if regeneration pseudo-cycles in the latter (see e.g. 
Ref. 1 34 1) share a common physical origin with nonlinear 
cycles in wall-bounded shear ffows 1 27J . 
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Appendix: Symmetries 



Nagata 1541 identified several possible symmetries for three-dimensional nonlinear hydrodynamic solutions in cen- 
trifugally unstable (Rayleigh-unstable) Taylor-Couette flow in the thin-gap limit, which corresponds to a cartesian plane 
Couette flow with walls, rotating along its span wise z axis. These symmetries can be adapted to the shearing box and 
generalized to IMHD flows. 

We consider a physical domain with < x < L^, Q < y < Ly and Q < z < L^. For visualization purposes, we assume 
that the observer is comoving with the base flow at x = Lxl2 and therefore introduce 



y =y- 



-fst. 



(A.1) 



In Fig. |4j the coordinates of the lower left comer of each box are x = 0, y' = 0, z = 0. The various fields entering the 
nonlinear dynamo cycle solution presented in this paper have the following generalized symmetry: 



Ur = 



Ux,ttif) cos (^^e^) sin (ky^y + kx(t)x^ -h Ux,oo(t) cos (kzoz) cos (kyoy + kx(t)x^ 
Uy^QQit) COS (^^e^) sin (ky^y -\- kx(t)x^ COS (^^o^) COS (kyoy -\- kx(t)x 

Uz,tt{t) sin (^^e^) COS {ky^y + kx(t)x^ + Uz,oo(t) sin (^^o^) sin [kyoy' + kx(t)x^ 



(A.2) 



B : 



Bx = 

By = B^ 



B: 



B. 



Bx,Qo(t) COS (kz^z) sin (kyoy + kx(t)x^ + 5x,oe(0 COS (^^o^) COS (ky^y -\- kx{t)x) 
By,Qo(t) COS (kz^z) sin (kyoy + kx(t)x^ + By^ot(t) cos (^^o^) cos (ky^y -h ^^(O-^) 
^,eo(0 sin (kzcz) COS (kyoy + kx(t)x^ + 5z,oe(0 sin (^^o^) sin (ky^y -h kx(t)x^ 



(A3) 



where the e and o subscripts indicate that the associated discrete wavenumbers are based on even and odd relative integers, 
respectively. kx(t) in each individual trigonometric expression is defined implicitly by Eq. Q using the ky wavenumber 
of the same expression. It can be checked that this symmetry is conserved by the IVIHD equations ([l)-([3' ^^en needed, 
we enforced numerically that the dynamical evolution took place in the invariant subspace defined by this symmetry by 
imposing it in the initial conditions and by further enforcing it every shearing time during the numerical integrations in 
order to avoid the growth of non-symmetric numerical noise. 



